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SUMMARY 

We provide a comparative study of the Subspace Projected Approximate Matrix method, abbreviated SPAM, 
which is a fairly recent iterative method to compute a few eigenvalues of a Hermitian matrix A. It falls in 
the category of inner-outer iteration methods and aims to save on the costs of matrix- vector products with 
A within its inner iteration. This is done by choosing an approximation Aq of A, and then, based on both 
A and Aq, to define a sequence {Aj.)]^^q of matrices that increasingly better approximate A as the process 
progresses. Then the matrix Aj^ is used in the kth inner iteration instead of A. 

In spite of its main idea being refreshingly new and interesting, SPAM has not yet been studied in detail by 
the numerical linear algebra community. We would like to change this by explaining the method, and to show 
that for certain special choices for Aq, SPAM turns out to be mathematically equivalent to known eigenvalue 
methods. More sophisticated approximations Aq turn SPAM into a boosted version of Lanczos, whereas it 
can also be interpreted as an attempt to enhance a certain instance of the preconditioned Jacobi-Davidson 
method. 

Numerical experiments are performed that are specifically tailored to illustrate certain aspects of SPAM 
and its variations. For experiments that test the practical performance of SPAM in comparison with other 
methods, we refer to other sources. The main conclusion is that SPAM provides a natural transition between 
the Lanczos method and one-step preconditioned Jacobi-Davidson. Copyright © 201 1 John Wiley & Sons, 
Ltd. 

Received . . . 

KEY WORDS: Hermitian eigenproblem; Ritz-Galerkin approximation; Subspace Projected Approxi- 
mate Matrix; Lanczos; Jacobi-Davidson. 



1. INTRODUCTION 

We provide a comparative study of SPAM [26]. SR\M, which stands for Subspace Projected 
Approximate Matrix, is a fairly recent (2001) method for the computation of eigenvalues of a large 
Hermitian matrix A. Like the Davidson method [8], SPAM was originally developed for matrices 
that arise from applications in Chemistry. It was only in [7], many years after its conception, that 
Davidson's method was given proper attention by the numerical linear algebra community. As far as 
we can tell, also SPAM has been neglected by the numerical linear algebra community. Moreover, 
even though a number of citations [6, 14, 16, 20, 33] within the Chemistry and Physics communities 
over the past years demonstrate awareness of its existence, no studies of its mathematical properties 
seems to exists. 

SPAM belongs to the category of inner-outer iteration methods and is interpreted by its inventors 
as a modification of the above mentioned Davidson method. It is based on the following observation. 
Even when sparse, the computational effort necessary to carry out a matrix-vector multiplication 
with A can be significant and often represents the bottleneck of the total computational effort. SPAM 
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reduces the costs of matrix-vector products with A by replacing its action within the inner iteration 
of the algorithm with a sparser or more structured approximation. By doing so, it attempts to slash 
the overall computational cost. The idea is not altogether new and is related to a certain type of 
preconditioning, called one-step approximation in the Jacobi-Davidson method. See Section 4.1 of 
[23]. There too, the matrix A is, in the inner iteration, replaced by a preconditioner. The originality 
of the approach in [26] lies in the fact that the action of the preconditioner is only applied to the 
subspace in which the action of A has not yet been computed in the outer iteration of the method. 
Consequently, the approximated action of A in the inner iterations is likely to become more and more 
accurate as the number of outer iterations increases. Note that nonetheless, only one approximate 
matrix Aq is needed. Intuitively, one would expect that SPAM would outperform Jacobi-Davidson 
with one-step approximation. 

Since the main idea of SPAM is refreshingly new and potentially interesting, with links to other 
eigenvalue methods and selection techniques, we would like to bring it to the attention of the 
numerical linear algebra community. Indeed, following an initial studies that led to a MSc thesis by 
the second author, there is a need to understand the method in more detail, not only in theory, but also 
to submit it to a more illustrative set of numerical experiments than those in [26] . In the experiments 
in [26], the overall efficiency of SPAM was the main interest instead of its mathematical position 
within the class of iterative eigensolvers. In particular, we would like to point out the similarities 
with and differences to strongly related and well-known iterative methods such as the Lanczos 
method [9] and the Jacobi-Davidson method of Sleijpen and van der Vorst [23]. 



1.1. Outline of the results in this paper 

We will show that for certain choices of the approximate matrix Aq, the SPAM method is 
mathematically equivalent to methods such as Lanczos [9] or the Riccati [4] method, another attempt 
to improve upon the Jacobi-Davidson [23] method. Further, we will see that a Schur complement- 
based choice for the action of A outside the Ritz-Galerkin subspace that is being built in the outer 
iteration naturally leads to a connection with harmonic Rayleigh-Ritz [3] methods. Next, we show 
that choosing Aq such that A — Aqis positive semi-definite has, at least in theory, an agreeable effect 
on the approximations obtained in the inner iteration in comparison to choosing Aq without such 
a restriction. Numerical experiments suggest that this also works through into the outer iteration. 
We comment on how such approximations Aq can be obtained in the context of discretized elliptic 
PDEs [5, 31] but also from a purely algebraic point of view. Finally, we present a variety of detailed 
numerical illustrations of the performance of the method in comparison with the Lanczos and the 
Jacobi-Davidson method. These illustrations do not merely aim to show that SPAM is a suitable 
method to solve eigenvalue problems (which was for a large part already taken care of in [26]), but to 
emphasize the role of approximations Aq from below and to show the similarities and discrepancies 
with Lanczos and Jacobi-Davidson, such that the SPAM method can be put into a proper perspective. 



2. SPAM AND SOME OTHER SUBSPACE METHODS 

Eigenvalue problems are among the most prolific topics in Numerical Linear Algebra [1, 13]. 
In particular, the continuous increase in matrix sizes turns the understanding of known iterative 
methods, as well as the development of more efficient ones, into an ample field of work within the 
global topic. Successful methods like the Implicitly Restarted Arnoldi [27] and the Krylov-Schur 
[28] methods and their symmetric counterparts are nowadays among the most competitive. For 
convenience and in order to set the ground for what is to come, we opt to outline the main concepts. 
For more detailed considerations on both theoretical and practical aspects of the numerical solution 
of large eigenvalue problems, we refer to [19, 12, 29, 30]. 
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2.1. Ritz values and vectors 

Throughout this paper, A is a Hermitian nx n matrix with eigenvalues 

An<An-i<---<A2<Ai, (1) 

Let F be an n X /c matrix with mutually orthonormal columns and define the /c-dimensional 
subspace V of by 

V = {Vy I yeC^}. (2) 

In the context of iterative methods, V is called the search subspace. Let be such that is 
unitary and write 

^ M i?* 



A = {V\V^rA{V\V^)= - - . (3) 
The eigenvalues of the k x k matrix M = V*AV, 

l^k < M/c-l < • • • < < Ml, (4) 

are called the Ritz values of A with respect to V. The vectors Ui = Vzi, where 2:1, . . . , is an 
orthonormal basis for consisting of eigenvectors of M belonging to the corresponding /i^, are 
the Ritz vectors of A in V. The residuals Vi = Aui — Uijii for the respective Ritz pairs (/i^, Ui) satisfy 

Aui - Uiiii = fi ± V. (5) 

Each Ritz pair {jii^ Ui) is also an eigenpair of the n x n rank-/c matrix VMV and is interpreted as 
an approximation of an eigenpair of A. See [15, 19, 21]. 



2.2. Rayleigh-Ritz and subspace expansion 

The Rayleigh-Ritz procedure is the first stage in iterative methods for eigenproblems and consists of 
computing the Ritz pairs from V. The computation of 5* in (3) is not needed, nor feasible for reasons 
of efficiency. However, a cheaply available by-product of the computation of AV is the matrix 
R = AV — VM, where R = VIR. Its columns are the respective residuals (5). In the second stage, 
the search subspace V is expanded. Different definitions of the expansion vector distinguish the 
different iterative methods. Each strategy results in a sequence of nested spaces 

Vl c V2 c • • • C Vn-l c Vn (6) 



and has the objective to obtain accurate Ritz pairs while spending only minimal computational 
effort. One of the strategies will result in SPAM. Other strategies lead to methods with which SPAM 
will be compared in this paper: the Lanczos [9], Jacobi-Davidson [23], and Riccati [4] methods. 



2.3. The Lanczos method 

The Lanczos method [9] defines Vj+i as Vj span{r} where r ± Vj is any of the current residuals 
from (5). This results in a well-defined method: starting with some initial vector vi with \\vi\\ = 1 
that spans the one-dimensional search space Vi , it can be easily verified by induction that regardless 
which residual is used for expansion, the sequence of search spaces that is defined, equals the 
sequence of Krylov subspaces 

Vj = JC^ {A,vi) = spd.n{vi, Avi, . . . , A^-^i}. (7) 

Due to (5), the matrix Vk with the property that its column span equals /C^(A, vi) can be chosen to 
have as columns vi and the normalized residuals 

^,+1 = ^, jG {!,..., fc}, (8) 
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where rj is a residual from the j-dimensional search space. From this the so-called Lanczos relation 
results, 

AVk = Vfc+iMfc,fc+i, with Mfc,fc+i ^ 



Pkel 



(9) 



Here, Mj. is a /c x A: tridiagonal matrix, e/e is the kth canonical basis vector of and pk is a scalar. 
The following trivial observation will have its counterpart in the discussion of the SPAM method in 
Section 2.4. 

Remark 2.1 

If the Lanczos method runs for the full number of n iterations, it produces a unitary matrix Vn 
that depends only on A and the start vector \\vi\\. The /cth leading principal submatrix of the n x n 
tridiagonal matrix M = V*AVn is the matrix Mk from (9) whose eigenvalues are the Ritz values 
after k iteration steps. 



2.4. The Sub space Projected Approximate Matrix ( SPAM) method 

In the Subspace Projected Approximate Matrix (SPAM) method [26], the expansion vector is a 
suitable eigenvector of an approximation of A. This approximation has a cheaper action than A 
itself. Thus, the matrix-vector products within the inner iteration that is needed to compute this 
eigenvector, will be cheaper than for instance in the Jacobi-Davidson [23] and Ricatti [4] methods. 
These methods, which we will explain in more detail in Section 2.5, both use A itself in their inner 
iteration. 

A central observation in [26] is that the action of A on V that has already been performed in 
the outer iteration, can be stored in a matrix W = AV, and be re-used within the inner iteration 
at relatively low costs. Thus, the action of A in the inner iteration only needs to be approximated 
partially. The resulting approximation is then different after each outer iteration step, even though 
only one approximate matrix Aq is provided. Its action is merely used on less and less of the total 
space. As such, SPAM may be interpreted as a discrete homotopy method. 



2.4.7. General description of SPAM. Let Aq be an approximation of A. In Section 3 we will 
conmient on how this approximation can be chosen. For now, assume that Aq is available and define, 
in view of (3) 

S = V^AoVl (10) 

and 

/Vf /?* 

Ak = {V\V^)Ak{V\V^r, where = ^ ^ . (11) 

The subscript k of A^ refers to the number of columns of V. In [26], the matrix A^ is called 
a subspace projected approximate matrix. This is motivated by the fact that A^V = AV and 
VAk = VA. In particular, since M = VAV = V*AkV, both Ak and A have the same Ritz pairs 
in V. This will be exploited to derive bounds for the eigenvalues of A^ in Section 3. This is of 
interest since the search space V in the outer iteration is expanded with an eigenvector of A^. Note 
that the action of Ak on does not equal the action of Aq on V^. Since A* = A, the action of Ak 
on equals, in fact, the action of A ink linearly independent functionals on C^. 
With the convention that Eq = /, write 

TVk = V^Vl =I-VV\ (12) 

This shows, together with equations (10) and (3), that 

Ak = -VV*AVV* + AW* + VV*A + TlkAQlik- (13) 

Thus, the action of Ak can benefit from the stored action of A on V as follows. With W — AV we 
have that 

AkV = -VMV*v + WV*v + VW*v + likAQlikV (14) 

where we have used (12) to avoid the numerically infeasible formation of V^. Note that if v G V^, 
the first two terms vanish. In view of Remark 2.1 we now observe the following. 
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Figure 1. Arrowhead updates from Ai^_i to 4/^ for consecutive values of k. 



Remark 2.2 

If SPAM runs for the full n iterations, it produces a unitary matrix Un that depends only on A 
and Aq. The kth leading principal submatrix of the n x n matrix M = U^AUn then contains the 
Rayleigh-Ritz approximations after k steps of the outer iteration. 

For theoretical purposes and without loss of generality, we assume that Vx in (11) contains 
precisely the basis for the orthogonal complement of V that SPAM is about to produce in future 
iterations. With respect to this basis, A^ is an update of A^-i of arrowhead type, in the sense that 





" 








— ^k-l = 


0* 


r 


r 




0* 


t 







(15) 



where each entry in the arrowhead formed by t G C^~^^r G M and t*, is the difference between the 
corresponding entries of A from (3) and Aq = (^| Vx)*Ao(F| V±) from (11). Thus, with respect to 
the basis defined by the columns of (F| Fx), the matrix Aq simply transforms step by step into A 
in the sense that after k steps, the first k columns and rows have changed into those of A, and the 
resulting matrix is called A^. This is visualized in Figure 1. 

On the original basis, this transformation is described in the next proposition. Note that in this 
proposition, v is the eigenvector of interest of orthonormalized toVk-i. 

Proposition 2.3 

Let k>l. Write {Vk-i\v\V±) = {V\V±), thus v = Vck. Then Ak is the following indefinite 
Hermitian rank- 2 update of Aj^-i, 



Ak = Ak-i + uv'' + v^x* = Ak-i + {u\v) 



1 

1 



where 



iu\vr, 



Uk-i - -vv* {A - Ao)v = [Uk + -vv* {A - Ao)v. 



I A , . . . /„ 1 



Proof. Combining (15) with iV\V±){A - Ao){V\V±)* = A-Aov/e find 

= (0 efe •■• Y {A-Ao)ek = {0\v\V^nA-Ao)v. 



el+ek{0* T t* )- eurel iVu-MV^)* 




(16) 



(17) 



(18) 



Therefore, substituting (18) into 

Ak - Ak-i = {Vu-i\v\V^) 




we arrive, using that v = {Vk-i \v\Vs_), at 

Ak - Ak-i = Ilk-i{A - Ao)vv* + vv*{A - Ao)nfc_i - vv*{A - Ao)vv*. (19) 
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Splitting the most right term in two equal parts and rearranging the terms proves the formula for u 
in (17). Since by (12) 

Uk-i-vv'' =Uk, (20) 
also the second equality is proved. □ 
Remark 2.4 

Note that the result of Proposition 2.3 can also be derived in an alternative way, starting from (13). 
The proof presented here also shows the validity of (15). 

From (3) and (1 1) we see that 

rank(A - Aj,) = rank(A - Aj,) < n - k, (21) 
and thus, even though — A^-i has rank-2, the update - or maybe better downdate 

A-Ak=A- Ak-i + {Ak-i - Ak) (22) 

will generically decrease the rank of A — A^-i with at most one. This remark goes hand in hand 
with the observation that even though the approximations A^ of A may seem unusual, the viewpoint 
of considering the reverse sequence = A — An^ A — A^-i, . . . ^ A — Ai^ A — Aq as increasingly 
better approximations of A — Aq is very natural indeed: they form a sequence of Rayleigh-Ritz 
approximations to A — Aq in the orthogonal complements of the spaces Vk. 

In the outer iteration, the products Av and V*Av were computed. Thus, both v*Av and U^^iAv 
are available in (17) without additional computational costs. Furthermore, since v is orthogonal to 
the {k — 1) -dimensional search space, we have that v*Ak-iv = v^Aqv. Now, because v is the result 
of orthogonalization of an eigenvector of A^-i to Vfc_i, also v'^Aqv can be retrieved from the inner 
iteration. Thus, the vectors u and v in the updating procedure (16) are cheaply available. 

Remark 2.5 

Of course, the update (16) itself should not be performed explicitly because it will generally result 
in fill-in of originally sparse matrices. 

2.4.2. Choice of the method for the inner iteration of SPAM. In [26], the authors suggest to use 
Davidson's method [8] to solve the eigenvalue problem for Ak, but of course any other method can 
be adopted. Apart from the Lanczos method, also the Generalized Davidson method from [17] was 
tested as inner method in [26]. Other possibilities include, for instance, the Jacobi-Davidson [23] 
method^ . The latter can be a good option because it often needs only a few iterations to converge 
if a good start vector is available. This start vector may be either the eigenvector approximation of 
Ak-i, or the current eigenvector approximation of A from the outer iteration. We will study this 
choice in Section 2.5. 

Remark 2.6 

SPAM should first of all perform well under the assumption that the eigenproblem for Ak is solved 
exactly. This will be investigated in the numerical illustrations in Section 4. 

In [26] it is also noted that SPAM itself can be chosen in the inner iteration. This leads to a 
recursive multilevel version of the method, and assumes that a whole sequence of approximating 
matrices of A is available, each having a cheaper action than its predecessor. The eigenvector 
computed at a given level is then used as expansion vector at the first higher level. 



^M. Hochstenbach presented this option at the 2006 GAMM/SIAM Apphed Linear Algebra conference. 
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2.5. Comparing SPAM with the Jacobi-Davidson and the Riccati method 

The philosophy of SPAM is to use all available (computed) information from (3), i.e., M, R and i?*, 
to determine the expansion vector, and consequently, that only S needs to be approximated. The 
Jacobi-Davidson [23] and Riccati [4] methods partially share this philosophy. Instead of R, they 
only use the residual corresponding to the selected eigenvector approximation. On the other hand, 
in their most simple forms, they do not approximate the matrix S but use its full action. In this 
section we will outline their similarities and differences. 

Remark 2. 7 

Since in the Lanczos method all residuals are linearly dependent, also the Lanczos method uses all 
information about the residual. However, Lanczos has no inner iteration. This will make sense from 
the point of view taken in Section 3. There we will show that for the choice Aq = 0, also SPAM 
needs no inner iteration, and that this choice makes SPAM mathematically equivalent to Lanczos. 

2.5.7. The Riccati and the Jacobi-Davidson methods. Given a Ritz pair (/i, u) with residual f ± V, 
generally each eigenvector of A has a multiple that equals u^t, with t ± a so-called orthogonal 
correction to u. Indeed, let X be such that {u\X) is unitary. Then with S = X*AX and f = Xr, 



A{u\X) = {u\X) 



fi r 
T S 



(23) 



and the orthogonal correction i equals Xp where p can be verified to satisfy the generalized 
algebraic Riccati equation 

(S — iJ^I)p = — r + pr*p. (24) 
Transforming (24) back to the original basis, shows that t solves 

t±u and {I -uu*){A- -uu*)t = -r^trt. (25) 

In [4], solutions of (26) were approximated by means of Rayleigh-Ritz projection in a ^-dimensional 
subspace U of u-^ and a suitable one was selected as expansion vector for V. This idea was intended 
as an enhancement of the Jacobi-Davidson method [23]. This method neglects the quadratic term 
tr*t from (25) and uses instead the unique solution i of 

i±u and {I-uu*){A-^I){I-uu*)i=-r (26) 

to expand V. Jacobi-Davidson is simpler than the Riccati method in the sense that only a linear 
system for t needs to be solved. It can be interpreted as an accelerated Newton method [25]. 
However, much more than the Riccati method, Jacobi-Davidson suffers from stagnation in case 
the term tf*t from (25) is not small. On the other hand, if tf*t is small enough, Jacobi-Davidson 
converges quadratically, as one would expect from a Newton type method. This shows that one 
should be careful in proposing alternatives to the correction equation (26). For instance, in [10], the 
authors investigated the effect of solving the following alternative correction equation, 

i±V and (/- VV*)(A-/i/)(/- VV*)t = -f. (27) 

At first sight this seems to make sense, because it directly looks for a correction orthogonal to 
V. Also, the conditioning of the linear equation (27) may be better than (26) in case the search 
space contains good approximations of eigenvectors belonging to eigenvalues close to /i. However, 
orthogonalizing t from (26) to V generally does not result in i from (27) and the price to pay is 
that p — t|| is not of higher order, as is ||t — Indeed, in [10] it is shown explicitly that, apart 
from some exceptional cases, the quadratic convergence of Jacobi-Davidson is lost, whereas in 
those exceptional cases, both expansions are equivalent. Numerical experiments in [10] confirm the 
above observations. 



Copyright © 2011 John Wiley & Sons, Ltd. 
Prepared using niaauth.cis 



Numer. Linear Algebra Appl. (2011) 
DOI: 10.1002/nla 



8 



J. H. BRANDTS AND R. REIS DA SILVA 



Remark 2.8 

The method using the correction equation (27) may at first sight also resemble SPAM. However, 
as we will see in the section to come, the use of V in SPAM is of a different nature than in (27), 
where, the correction is sought orthogonal to V. In SPAM, it is sought in the whole space and only 
orthogonalized to V afterwards. 



2.5.2. Preconditioning in the Jacobi-Davidson method. In the original paper [23] on the Jacobi- 
Davidson method, reprinted as [24], preconditioning is discussed as follows. Suppose that an 
approximation of A is available. It is shown in [23] how to apply such a preconditioner to (26), 
which is a linear equation, though not with system matrix A. To be explicit, since (/ — uu*)i = £, 
we have that 

{A - = -eu - r. (28) 
where e is such that t ± u. Or equivalently, written as an augmented system. 



(29) 



u 




" t ' 




—r 







e 








Thus, an approximation fo of U together with an approximation sq of e can be obtained simply by 
replacing A by Aq in (29). The pair fo, can be computed as 



to = -Sq^Aq - III) - (Ao - III) ^f, with So = 



^i*(Ao-/i/)-^f 
u*{Aq - ilI)-^u' 



(30) 



This approximation is called a one step approximation in [23]. It was observed that setting £0 = 
in (30), the Davidson method [8] results. With Aq = A, which corresponds to Jacobi-Davidson with 
full accuracy solution of the correction equation, (30) becomes 



t 



-e{A — III) — 



(31) 



and since i is then orthogonalized to u, the method is mathematically equivalent to an accelerated 
shift and invert iteration that works with {A — iil)~'^u. It is argued, and demonstrated by 
experiments in [23], that Jacobi-Davidson combines the best of those two methods. Of course, a 
natural next stage in preconditioning is to use the matrix 



AU _ 
^0 — 



Aq — III u 





(32) 



as preconditioner within the iterative method that aims to solve (29). In each step of such a method 
one would need to solve a system with Aq. This can be done by solving two systems as in (29)- 
(30) in the first step of the inner iteration. In each consecutive step, only one system of the form 
(Ao — iil)z = y would need to be solved. 



2.5.3. One step approximation of the SPAM eigenproblem for A^. In the SPAM method, the 
expansion vector for the Ritz Galerkin subspace in the outer iteration is a relevant eigenvector 
Vk of A^. In principle, any eigenvalue method can be used to compute an approximation for v^, 
but observe that the starting point is as follows. In the outer iteration, we have just solved a /c x /c 
eigenproblem for M = V'^AV, and a Ritz pair (/i, u) with \\u\\ = 1 and with residual r = Au — iiu 
has been selected. The matrix A^ is now available, either explicitly or implicitly. Since A^V = AV 
and V*Ak = V*A, the Ritz pair (/i, u) for A with respect to the current search space Vk is also a 
Ritz pair for A^. Thus we can exploit the fact that Vk contains, by definition, good approximations 
of the relevant eigenvectors of Aj with j < k, and use it as initial search space for a Ritz Galerkin 
method applied to Ak to approximate Since Vk is generally not a Krylov subspace, the Lanczos 
method is not a feasible candidate. The Jacobi-Davidson method is. The correction equation for the 
first step of the Jacobi-Davidson method in the inner iteration can be set up without any additional 
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computations: 



Ak- III u 






tk 




—r 




_ £k _ 








(33) 



Since quadratic convergence in the outer iteration cannot be expected even if an exact eigenvector 
of Ak would be computed, we study the effect of applying only one iteration of Jacobi-Davidson 
in the inner iteration. This is also motivated by the fact that the initial search space Vk for Jacobi- 
Davidson applied io Ak may be relatively good and may result in quadratic convergence in the inner 
iteration. 

Remark 2.9 

If only one step of Jacobi-Davidson is applied, then after solving from (33), the new 
approximation v for the eigenvector Vk of Ak would lie in the space Vk {tk)- It would not be 
necessary to actually compute this approximation, because 



Vk(B{v) <zVk(B{tk). 



(34) 



Thus, instead of computing the eigendata ofa(/c + l)x(/c + l) matrix, the expansion of the Ritz- 
Galerkin space in the outer iteration can be done immediately with tk. 

SPAM, in which the eigenproblem for A/e is approximated with one iteration of Jacobi-Davidson, 
we will refer to as one step SPAM, abbreviated by SPAM(l). We will talk about Full SPAM if the 
eigenproblem for is solved to full precision. 

2.5.4. Comparing one step Jacobi-Davidson with SPAM(l). SPAM(l) can best be compared with 
preconditioned Jacobi-Davidson with one step approximation, as described in Section 2.5.2. The 
only difference between the two method is that in iteration k of SPAM(l) the preconditioner Ak is 
used, whereas in one-step Jacobi-Davidson this is Aq. As such, SPAM(l) can be seen as an attempt 
to enhance this type of preconditioned Jacobi-Davidson. We will now investigate the effect of this 
attempt. 

Lemma 2.10 

Assume that Vk+i = V/c {v) with v ^Vk and \\v 



1. Then 

Ak+i = -vv'^Avv'^ + Aw'' + vv'^A + (/ - vv'')Ak{I - v'^*). 
Proof. By substitution of the defining relation for Ak. 



(35) 
□ 



Corollary 2.11 

Let Vi be the span of the relevant eigenvector u of Aq and let ji = u'^Au and f = Au — jiu. Then 
the solution ti from the system 



(36) 



in the first iteration of SPAM(l), to be used to expand Vi, coincides with the solution to from the 
system 

Aq — jll u 

u* 



Ai — jil u 








—r 


ii* 














" to ' 




— f 




_ £o _ 








(37) 



solved in the first iteration of Jacobi-Davidson with one-step approximation using Aq. 
Proof. The linear system (33) for SPAM(l) with A: = 1 is equivalent to 

ti ± {I — uu'^){Ai — /il){l — uu*)ti = — f, 



(38) 



where u is 3. unit vector spanning Vi. Substituting the expression (35) for Ai immediately proves 
the statement. □ 
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Thus, there is no difference in the first iteration. There is, however, a difference in further 
iterations. To study this difference we take a different viewpoint. Above, we approximated the 
SPAM inner eigenvalue problem by a linear correction equation which made it suitable for 
comparison with one step Jacobi-Davidson. The opposite viewpoint is also possible, which is 
to interpret the Jacobi-Davidson correction equation with one- step approximation as an exact 
correction equation of a perturbed eigenproblem. 

Lemma 2.12 

Given u with \\u\\ = 1 and /i = u'^Au and r = Au — jiu. Define for a given approximation Aq of A 
the matrix 

Au := -uu'Auu' + uu'A + Auu" + (/ - uu')Aq{I - uu*) 

= uiiu" + ixr* + TO* + (/ - uu"')Aq{I - uu"). (39) 

Then (^/, /i) is a Ritz pair of Au in the one-dimensional span of u with residual r = AuU — jiu and 
with the equation 

(/ - uu*){Aq - - uu*)t = -r (40) 

as its exact (i.e., without preconditioning) Jacobi-Davidson correction equation. 

Proof. It is easily verified that u*AuU = fi and AuU — fiu = r and thus (/i, u) is a Ritz pair for Au 
with residual r. Since moreover, 

(/ - uu'')Au{I - uu'') = (/ - ixix*)Ao(/ - uu*) (41) 

its correction equation is precisely (40), or equivalently, (30). □ 

Note that Au from (39) is the subspace projected approximate matrix for the one dimensional 
subspace spanned by u. Now, with 

Vk=U®{u), (42) 

where U is the orthogonal complement of the span of the relevant Ritz vector u in the current search 
space Vfc, we have, similar as in (35) and (39), that 

Ak = uiiu* + ur* + ru* + (/ - uu*)Au{I - uu*). (43) 

Here, Au is the subspace projected approximated matrix corresponding to the subspace U. Now, 
the opposite viewpoint mentioned above, is the observation that in one step Jacobi-Davidson, the 
expansion vector is (an arguably good approximation of) the relevant eigenvector of Au in (39), 
whereas in SR\M, it is (an arguably good approximation of) the relevant eigenvector of A/e in (43). 
Both matrices have now an appearance that is suitable for studying their differences and similarities. 

The most important observation is that neither correction will lead to the unique correction that 
results in quadratic convergence in the outer iteration. Second, since both matrices Au and A^ differ 
only in their restriction to the orthogonal complement of u, the difference in the methods they 
represent will be marginal if the residual is already small. Since, as already mentioned in Corollary 
2.11, not only at the start but also in the first iteration both methods coincide, the difference from 
the second iteration onwards will probably be very small, especially if Aq provides a good initial 
approximation of the relevant eigenvector. Finally, even though Au uses less information of A than 
Ak, it does use the optimal information in some sense. It may even be a disadvantage to use more 
information, because this involves approximations of eigenvectors orthogonal to the eigenvector 
of interest. The above observations will be tested, and confirmed by our numerical illustrations in 
Section 4. 



3. SELECTING THE MATRIX Aq IN THE SR\M METHOD 

Here we study some of the effects of the choice for Aq on the iterative approximation process 
in SPAM. We will assume for simplicity that the largest eigenvalue of A is the target, although 
replacing Ahy - A, we might as well have set the smallest eigenvalue of A as a target. 
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3.1. Some simple choices for Aq 

It is instructive to study the consequences of the choice = 0. This generic parameter- free choice 
may seem dubious at first sight, but it is not. First note that since the start vector of SPAM is the 
relevant eigenvector of = 0, this is necessarily a random vector v e C^, with \\v\\ = 1. Thus, we 
set Vi = spanji;}. Write /i = v*Av and r = Av — vji. Then, with such that (v| V±) is orthogonal, 
and r = y^r we find that 



A = {v\V^) 



r 



S 



(44) 



and consequently, replacing S by the zero matrix, the next approximating matrix Ai from (11) is 
defined by 



= {v\V^) 



r 



r 




(45) 



As shown already in Proposition 2.3, A i is a simple rank- two matrix, that on the basis defined by 
the columns of (v| Vx) is of arrow-head type. It has two nontrivial eigenpairs Aiw± = 9±w±, where 



w± 



e^v^r and ^± = ^/i ± W ^/i^ + ||r||2. 



(46) 



Assuming that A is positive definite will lead to the selection of Wj^ for the expansion of the search 
space. Since Wj^ is a linear combinations of v and r, we find that 



V2 = spanj'u, r} = K?{A, v), 



(47) 



and the two eigenvalue approximations computed in the outer loop of the SPAM method are the 
same as in the Lanczos method. This is, of course, not a coincidence. 

Theorem 3.1 

If the goal is to find the largest eigenvalue of a positive definite matrix A, the SPAM method with 
= is mathematically equivalent to the Lanczos method. 

Proof. Let Aq = 0. Then the eigenvalues of A ^ in (1 1) are those of the n x n Hermitian arrowhead 

(48) 



M R* 
R 



The Cauchy Interlace Theorem immediately gives that the largest k of them are each larger than or 
equal to the corresponding eigenvalue of M. This assures that the eigenvector of Ak that is selected 
for expansion is not from its null space but from its column span. From (13) we see that if = 
this column span is the span of V and AV. A simple induction argument shows that this span equals 

Remark 3.2 

If A is indefinite and we would like to find the eigenvalue closest to zero, the choice Aq = would 
lead to expansion vectors from the null space of Aq, and the method would be worthless. As a 
solution, A may be shifted by a suitable multiple a times the identity / to make it positive semi- 
definite. Equivalently, instead of using = we may choose Aq = al as approximating matrix. 
In both cases, it is easy to verify that SPAM will still be equal to Lanczos. 

The observant reader may have noticed that a peculiar situation has arisen. In the inner loop of 
SPAM, a better approximation of the largest eigenvalue of A was computed than the Ritz values 
from the outer loop. In view of the philosophy of inner-outer iterations, this in itself is not out of 
the ordinary, but its computation did not require any additional matrix- vector multiplication with A, 
nor with an elaborate approximation Aq of A. The following proposition, which uses the notation 
(1) and (4), makes this explicit. 
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Proposition 3.3 

With 6>+ as in (46) we have that fi < \\Av\\ < 0^. Moreover, if A is positive semi-definite, we find 
that 6>+ < A+, where A+ is the largest eigenvalue of A. 



Proof. Because r ± v and Av 



IJLV 
l|2 



r, by Pythagoras' Theorem we have 



.^'^\\Tf = \\Avf 



hence ^ < \\Av\\. Squaring 9+ from (46) gives 

1 



WAV 



- Li ^ jl 



(49) 



(50) 



which together with (49) shows that \\Av\\ < 6>+. Since S in (44) is positive definite whenever A is, 
combining (44) and (45) with Weyl's bound yields 6>+ < A+. □ 

The key issue here is that the inner eigenvalue approximations are much related to the so-called 
harmonic Ritz values [3, 21] of A. Indeed, assuming that A itself is positive semi-definite, these are 
the k positive eigenvalues of the at most rank-2/c matrix 



M 

R T 



where T = RM-^R\ 



(51) 



They can be computed without additional matrix vector multiplications with A. Note that harmonic 
Ritz values are usually introduced as the reciprocals of the Rayleigh-Ritz approximations of A~^ in 
the space AV. It is well known that for positive semi-definite matrices A, the harmonic Ritz values 
are better approximations of the larger eigenvalues of A than the standard Ritz values. We provide 
the short argument in Lemma 3.5. See also [3]. 



Proposition 3.4 

The matrix A^ can be decomposed as 



The blocks 



■ M 


R* ' 




R 


T 






' M ' 






R 


a 



M 
R 



M-^ [M R* ]. 



-M-^R* 
I 



(52) 



(53) 



span the range and null space, respectively, and the nonzero eigenvalues are the eigenvalues of the 
k X k matrix 

" M 



UM-^U* where 



R 



= QU 



(54) 



is a QR-decomposition. In particular, those eigenvalues are positive. 



Proof. The statements are all easy to verify. The positivity of the k nonzero eigenvalues follows 
from Sylvester's Theorem of Inertia. □ 

The k eigenpairs (Oj.Wj) of Aj. in (51) with positive eigenvalues we label as 



0<^fc<^fc_l<---<^2<^l. 



(55) 



The proposition shows that they can be computed by solving di k x k eigenproblem and that no 
additional matrix- vector products with A are needed. We can now easily prove the following bounds. 
See also [3, 19]. 



Lemma 3.5 

For all j G {1, . . . , /c}, we have that 



(56) 
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Proof. The left inequalities follow from the Cauchy Interlace Theorem applied io A^. Now, with 
A as in (3), recognizing the Schur complement A/M shows that, 



A -A. 



0* 
A/M 



(57) 



is positive semi-definite, hence the right inequalities follow from Weyl's bound. □ 

Observe that, as was the case with the choice Aq = 0, assuming that from the start of the SPAM 
method the eigenvector wi belonging to Oi was selected for expansion, the eigenvectors wj, called 
the harmonic Ritz vectors, lie in the column span of AV and hence, in /C^+^(A, v). 

Thus, even though we may have improved the inner loop eigenvalue approximations, the SPAM 
method is still equal to the Lanczos method. It does give, however, valuable insight into SPAM: 
Lanczos results from the choice Aq = 0, even after modification of Ak into the positive semi-definite 
matrix Aj^. In order to get a method different than Lanczos, we should use less trivial approximations 
that are based on the structure and properties of A itself, while aiming to retain similar inequalities 
as in Lemma 3.5. For this, we would need the matrices A — A^io be positive semi-definite. We will 
investigate this in the following sections. 



3.2. One-sided approximations 

Having seen that the trivial choice = 0, even after after a correction that turns all approximate 
matrices positive semi-definite, will generally lead to the Lanczos method, we now turn our attention 
to approximations from below, by which we mean Aq such that A — AqIS positive semi-definite. 

Lemma 3.6 

If Aq approximates A from below, then so does each matrix A^. 
Proof. Combining (3) and (10) with (1 1) we see that for all x G C^, 

- Ak)x X* 

where the last inequality holds because A — Aq is positive semi-definite. And thus. A — Aj^ is 
positive semi-definite, and hence, so is A — Ak. □ 

By Proposition 2.3, A^ — A^-i is an indefinite rank-2 matrix, hence it will generally not be true 
that Ak-i approximates Ak from below. 

Lemma 3. 7 

The following inequalities are valid generally, 

Oj^i < iij < Oj for all j G {1, . . . , k}, (59) 

whereas, if Aq approximates A from below, additionally 

Oj < Xj for all j G {1, . . . , n}. (60) 

Proof. The first set of inequalities applies because Ak has the same Ritz values as A. See also 
Section 2.4. 1 . It is well known that the Ritz values interlace the exact eigenvalues [19]. Since A — Ak 
is positive semi-definite due to Lemma 3.6, the equality Ak -\- {A — Ak) = A together with Weyl's 
bound [19] proves the second set of inequalities. □. 

Lemma 3.7 shows that if Aq approximates A from below, the approximations for the larger 
eigenvalues of A that are produced in the inner iteration, will never be worse than the ones obtained 
in the outer iteration. Moreover, they will never be larger than the corresponding exact eigenvalues. 
Thus, it indeed makes sense to expand the search space with the eigenvector that is computed in the 
inner iteration. Question that remains is how to obtain matrices Aq that approximate A from below. 



0* 

s-s 



x = x''Vl{A-Ao)V^x>0 



(58) 
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. = 3645 nz = 968 



Figure 2. Approximating a matrix A from structural engineering from below by a sparser matrix Aq by 
subtracting from A a definite matrix H; the sparsity plot of A (with 3648 nonzero entries) is on the left and 
of Aq (with 968 nonzero entries) on the right. See the experiments with this pair of matrices in Section 4.4. 



3.2.1. Algebraic construction of approximations from below. Clearly, for any positive definite 
matrix H we have that Aq = A — H approximates A from below, even though Aq itself may not 
be positive definite. The problem is of course how to choose H such that Aq is close to A in an 
appropriate sense, while its action is considerably less expensive than that of A. 

If A itself is a positive definite matrix, a purely algebraic option is at hand. Given an index set 
1 C {1, . . . , n} of cardinality m, let Ex be the matrix with the standard basis vectors e^, z G X as 
columns. Then set 

Ao=A-H, where H = E^E^AE^E^. (61) 

The matrix H is the result of a Rayleigh-Ritz procedure with as search space the column span of 
Ex. For a randomly chosen index set, this space has no a priori relation with A and thus H is 
probably a relatively poor approximation of A in comparison with, for example, a Krylov subspace 
approximation. 

Remark 3.8 

In this particular situation it is an advantage if H does not approximate A very well, because Aq 
should be close to A, not H. Notice also that Aq has zero entries at positions (z, j) for all i^j G X, 
and is thus always more sparse than A. A priori knowledge of A may lead to a more sophisticated 
choice of the index set(s) X. If the goal is to approximate the largest eigenvalues of A, the index set X 
could be chosen such that the smallest diagonal entries of A are selected to put in H. Consequently, 
Aq will share with A the largest diagonal entries, and this may increase its approximation quality. 
This is illustrated in Figure 2. 

Remark 3.9 

Notice that rank(74o) < 2{n — m). Especially for large m this may greatly simplify the computation 
of the eigendata of Aq and of Aj^ for small values of k. 

Remark 3.10 

Although A itself is positive definite, generally is not. It can be made positive definite by 
adding a Schur complement at the position of the zero block, as was done in (51). Since the Schur 
complement block is the smallest correction to Aq that makes it positive semi-definite, the result 
would still approximate A from below. However, its computation involves the evaluation of the 
inverse of an {n — m) x {n — m) matrix. It is not clear if the additional computational effort is well 
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spent. Also without the Schur complement, Aq has as many positive as negative eigenvalues and is 
likely to perform better than Aq = 0. 

The rather crude approach that we just described is easy to implement, and it can be applied in 
the context of the multilevel version of SPAM mentioned in Section 2.4.2: the solution of the inner 
iteration is done using SPAM itself, with an approximation of Aq that is based on a larger index set 
than the one that was used to construct Aq itself. 

3.2.2. Natural construction of approximations from below. A situation in which approximations 
from below are naturally available is the setting of discretized partial differential equations including 
an elliptic term, either by the finite difference method or the finite element method [5]. Then 
removing the positive definite discrete elliptic operator, either completely or partially, from the 
total discrete operator, results in an approximation from below. Indeed, let 1] c be a domain, 
and consider as an example the problem of finding eigenmodes for the linear operator L, defined by 



L{u) = Xu where L{u) = —ediv{JCVu) + cu, and u = on dQ. 



(62) 



Here e > is a parameter, JC : ft ^ M'^^'^(IR) maps into the symmetric positive definite matrices, 
and c e C{ft) is nonnegative. Discretizing this equation with the finite difference method will lead 
to an algebraic eigenvalue problem of the form sKx + Mx = ^x. The matrix K that represents the 
discretized diffusion is positive definite. Although sparse, it will generally have more fill-in than the 
diagonal matrix M that represents the reaction term, and, if e is small enough, have smaller entries. 
Thus, the total discretized operator A = K -\- M has = M as a candidate for the approximation 
matrix: its action is cheaper than the action of A and A — Aq = K is positive definite. A similar 
strategy can also be employed in the finite element method, when so-called mass lumping is used in 
the assembly of the mass matrices. 

Remark 3.11 

In this context, the algebraic method can be used to approximate the smallest eigenvalue of A by 
applying it to al — A with a such that al — A is positive semi-definite. Even though the largest 
eigenvalue usually has no physical relevance, together they would provide good estimates for the 
condition number of the A, which is indeed of interest. 



3.3. Cutting off the bandwidth 

Here we describe an obvious choice of approximating matrices that was used in [26] to illustrate 
the effectiveness of their method. It concerns symmetric banded matrices. Apart from the approach 
that we will now describe, in the numerical illustrations of Section 4 we also intend to apply our 
algebraic approximations from below to these matrices. 

Given < e < 1 and I < q < n — 1, define a matrix A = (Aij) by 



i ifi= j\ 

ifl<|i-j|<^; 
otherwise. 



(63) 



In [26] it is proposed to choose as approximating matrix for A a matrix A^ of the same type with a 
smaller half-bandwidth q < q. For instance. 



A = 



e 
2 

e 

^2 



e 
e 
3 
e 



then with Aq 



For each eigenvalue 6 of Aq there is an eigenvalue A of A with 



< 



/n. 



(64) 
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Indeed, the difference A - AqIS zero except for the bands qo ^lio q. Each non-zero row contains 
at most the numbers e^°+^ to e^, and (64) follows from the the Bauer-Fike [2] theorem and the finite 
geometric series sum formula. Thus, for values of e small enough, Aq may be a good candidate 
to approximate A even though it is generally not an approximation from below. Nevertheless, its 
eigenvalues are close to the eigenvalues of A and its action is cheaper than the action of A. The 
number of floating point operations required to compute Av is approximately twice as much as 
for an approximation Aq having half the bandwidth of A. In other words, each two matrix- vector 
products Aq are approximately equally costly as a single product with A. 

Cutting off the bandwidth in this fashion makes sense especially if the decay of the size of the 
off-diagonal entries in relation to their distance to the main diagonal is quick enough. Apart from 
the example above, this is the case in many applications where the boundary element method [31] 
is used to approximate integral equations. Notice that for such applications also the approach from 
Section 3.2.1 can be applied, resulting in both sparse and low rank approximating matrices Aq. 

4. NUMERICAL ILLUSTRATIONS 

In the previous sections, we have compared SPAM with both the Jacobi-Davidson method and the 
Lanczos method. We have also studied ways to define an appropriate approximate matrix Aq in the 
context of the SR\M method. Since the choice Aq = will effectively lead to the Lanczos method, 
our starting point in the upcoming numerical illustrations will be to consider SPAM as a boosted 
version of the Lanczos method. This is, of course, particularly justified if Aq is an approximation of 
A from below. As discussed in Section 2.5, SPAM can also be considered as an attempt to enhance 
the Jacobi-Davidson method with one-step approximation as preconditioning. Therefore, we will 
also present a comparison of SPAM and this version of Jacobi-Davidson. We end this section with 
a discussion of the numerical results. 

4.1. Objectives 

First, we list the methods and abbreviations that we use to describe them. 

• Lanczos (see Section 2.3); 

• JD(^): Jacobi-Davidson using £ steps of MinRES to approximate the solution of the exact 
correction equation (29) in augmented form; 

• JD(1,£): Jacobi-Davidson with one step approximation as preconditioner (see Section 2.5.2) 
using the matrix Aq, and £ steps of MinRES to approximate the solution of the correction 
equation (29), with A replaced by Aq; 

• Full SPAM: each eigenproblem for A^ solved in full precision; 

• SPAM(l): eigenproblem for Aj. approximated with one step of the Jacobi-Davidson method, 
correction equation (33) in augmented form solved to full precision (see Section 2.5.3); 

• SPAM(1,^): using £ steps of MinRES to approximate the solution of the correction equation 
for SPAM(l). 

Remark 4.1 

To minimize the size of the legends in the pictures, we sometimes write LZS for Lanczos, ESP for 
Full SPAM, SP(1) for SPAM(l), S12 for SPAM(1,2), JD13 for JD(1,3), etcetera. 

In the experiments, we will give illustrations of the following aspects of SPAM. 

• When a nonzero approximation Aq of A from below is used, less outer iterations of Full SPAM 
are needed to arrive close to the dominant eigenvalue of A than with the choice = 0, which 
is equivalent to the Lanczos method with a random start vector. 

• Even if the Lanczos method is started with the same approximation of the dominant 
eigenvector of Aq as Full SPAM, Full SPAM method will still outperform Lanczos in terms 
of the number of outer iterations. 



Copyright © 2011 John Wiley & Sons, Ltd. 
Prepared using niaauth.cis 



Numer. Linear Algebra Appl. (2011) 
DOI: 10.1002/nla 



ON THE SPAM METHOD 



17 



• Also for other eigenvalues, Full SPAM outperforms Lanczos; this may be expected because 
in Full SPAM the Ritz Galerkin subspace will be forced in the direction of the appropriate 
eigenvector in every iteration. In Lanczos, this is done only by means of the start vector. 
Of course, Lanczos allows efficient implicit restart strategies [27], but also Full SPAM may 
be restarted. We feel that the comparison would become diffuse and thus we refrain from 
incorporating restarts. 

• We investigate the effect of approximating the desired eigenvector of Ak with just one step of 
the Jacobi-Davidson method, i.e., we will be comparing Full SPAM with SPAM(l). 

• SPAM(1,^) will be compared with JD(1,^), both started with the same initial vector, i.e., the 
relevant eigenvector of Aq; i.e., both methods will spend the same number £ of matrix- vector 
products in their inner iteration, where in SPAM(1,^), the matrix will be A^, and in JD(1,^), 
the matrix will be Aq . From the viewpoint of this paper, this is the comparison that is of most 
interest. Not only is JD(1,^) the closest related to SPAM(1,^), the difference between the two 
is solely the fact that the action of A from the outer loop is taken into the inner iteration of 
SPAM(1,^), whereas in JD(1,^), this is not done. See the discussion in, and particularly at the 
end of Section 2.5.4. 

• Finally, we compare SPAM(1,^) with JD(^). This is perhaps the comparison that the authors of 
SPAM [26] had in mind: in the inner iteration of JD(£) the original matrix A is used, whereas 
in SPAM(1,^) it will be Aj,. 

We will comment on the computational costs of an inner iteration in comparison to having no 
such costs (as in Lanczos) or the full costs (Jacobi-Davidson), although these costs may depend very 
much on the specific problem and the available approximations and even on hardware parameters 
like available memory. 

4.2. Lanczos versus full SPAM: Reaction-Diffusion problem, various eigenvalues 

In this section we will compare Lanczos with Full SPAM. Our comparison is, for the time being, 
only in terms of the number of outer iterations. A first naive comparison uses a random start vector 
for Lanczos, but from then onwards, we will start Lanczos with the appropriate eigenvector of A^. 
The approximate matrix Aq will be contructed using both approaches described in Sections 3.2.1 
and 3.2.2. For this, we discretized the one-dimensional version of (62) on 1] = [0, 1] using finite 
differences with grid size h = 1/33 with the parameters e = ^,JC = 1, and c{x) = x{l — x)e^^. 
The resulting 32 x 32 algebraic eigenproblem is of the form Ax = Ax, where A = D-\-Ris the sum 
of the tridiagonal discretized diffusion D and the diagonal discretized reaction R. With the approach 
from Section 3.2.2 we approximate the largest eigenvalue and with the approach from Section 3.2.1 
the smallest eigenvalue. 

Natural approximation from below The left picture of Figure 3 illustrates the typical 
convergence of the Lanczos method with a random start vector. We display the k Ritz values at outer 
iteration k as circles above the value k of the iteration number on the horizontal axis. Due to the 
interlacing property, eigenvalue approximations "converge" to either side of the spectrum, and this 
is emphasized by the connecting lines between Ritz values for different values of /c, both upwards 
and downwards. In view of Theorem 3.1 we could also say that this picture belongs to SPAM with 
choice Aq = 0. In the right picture of Figure 3, we show in a similar fashion the convergence of 
SPAM, using Aq = R sls the approximate matrix, as suggested in Section 3.2.2. We see that the 
convergence towards the largest eigenvalues is stimulated. The costs for this faster convergence is 
solving a diagonal plus rank 2k eigenproblem in iteration step k. There exist efficient methods for 
such eigenproblems based on the secular equation and Newton's method. 

Algebraic approximation from below We also tested the algebraic approach from Section 3.2.1 
to construct approximations from below. We created a rank- 12 approximation Aq of A based on the 
largest diagonal elements of A. Full SPAM and Lanczos were started with the dominant eigenvector 
of Aq in order to approximate its dominant eigenvalue. The leftmost picture in Figure 4 shows that 
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the incorporation of an approximation Aq into Full SPAM has an effect that carries beyond that of 
creating a better start vector for Lanczos. In the rightmost picture of Figure 4, the same was done 
for the matrix 6/ — A, which is positive definite and has the smallest eigenvalue of A as dominant 
eigenvalue. Also here, SPAM outperforms Lanczos in terms of the number of outer iterations. In 
the middle two pictures of Figure 4, we plotted the absolute error in the second and fifth largest 
eigenvalue of A. The results are shown from iteration 2 and 5 onwards, respectively, because the 
Ritz-Galerkin method produces approximations of the second and fifth largest eigenvalues from that 
iteration onwards. Again, Full SPAM clearly outperforms Lanczos in both cases. Note that, when 
using Full SPAM to approximate the p-th largest eigenvalue, the Ritz-Galerkin subspace in the outer 
iteration is in each step expanded with the eigenvector belonging to the p-th largest eigenvalue of 
Ak. For a fair comparison, we also started Lanczos with the eigenvector of Aq belonging to the p-ih 
largest eigenvalue. 



LANCZOS METHOD WITH RANDOM START VECTOR SPAM WITH THE REACTION TERM AS APPROXIMATING MATRIX AD 
6| : . ^ ^ . . 1 6| ^ , ^ , , ^ 




qI 1 1 i ^ 1 1 1 gl ^ ^ ^ ^ ^ L_ 

5 10 15 20 25 30 5 10 15 20 25 30 



Figure 3. Lanczos method (left) with a random start vector, versus SPAM (right) with the discretized reaction 
term as approximating matrix Aq, and the largest eigenvalue of A as target. 




Figure 4. From left to right: Lanczos and Full SPAM approximating the largest, the second, the fifth, and the 
smallest eigenvalue of A, using algebraic rank- 12 approximation from below (for the smallest eigenvalue, 
we applied both the methods to 61 — A). Lanczos and Full SPAM used in each experiment the same start 

vector. 
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4.3. Lanczos versus Full SPAM and SPAM( 1 ): Banded matrices 

In this section we not only compare Lanczos with Full SPAM, but also with SPAM(l), by which 
we mean that the eigenproblem for in Full SPAM is approximated using one step of the Jacobi- 
Davidson method, as explained in Section 2.5.3. The correction equation that results is still solved to 
full precision. For our of experiments, we took the banded matrix from Section 3.3 of size 32 x 32 
and with q = 5 and e = 0.5. In [26], it was suggested to take for Aq the matrix A with a number of 
its outer diagonals put to zero. We repeated that experiment with Aq equal to the diagonal of A, such 
that here too, A^ is diagonal plus a rank-2/c perturbation. The comparison between Lanczos, Full 
SPAM and SPAM(l) is depicted in the left graph in Figure 5. This comparison is initially in favor 
of Lanczos. This may be due to the fact that the difference A — Aq is indefinite: it has only eleven 
nonnegative eigenvalues. In the middle left picture we took the tridiagonal part as approximation. 
This approximation is indeed positive definite and gives better results. Taking the simple algebraic 
approximation Aq from below from Section 3.2.1 based on the three largest diagonal entries of A, 
which is of rank 6, gives comparable results, but its low rank and higher sparsity make this choice 
more interesting. In the right graph in Figure 5, the largest eigenvalue of al — A was approximated 
with the approximation from below that kept the largest three diagonal entries of al — A, and 
thus the smallest three diagonal entries of A. In all cases, the positive effect of incorporating an 
approximation Aq goes beyond delivering a good start vector for Lanczos. Also in all cases, there is 
virtually no difference between Full SPAM and SPAM(l). 




Figure 5. Lanczos versus Full SPAM and SPAM(l) with diagonal (left), tridiagonal (middle left), and with 
algebraic rank-6 approximation from below (both middle right and right). In the three leftmost pictures, the 
target was the largest eigenvalue, at the right it was the smallest eigenvalue (i.e., the largest of al — A). 



4.4. Lanczos versus Full SPAM and SPAM( 1 ): matrices from structural engineering 

As a next set of experiments, we took some matrices from the Harwell-Boeing collection of test 
matrices. They have their origin in the area of structural engineering and are called bcsstk04, 
bcsstkO? and bcsstklO and have respective sizes 132 x 132,420 x 420 and 1024 x 1024. As 
approximating matrix Aq we took approximations from below keeping respectively the largest 
12, 20 and 180 diagonal entries. Recall that in Figure 2 we displayed the sparsity plots of A and 
Aq for bcsstk04.mtx. As was the case for the banded matrix in the previous section. Full SPAM 
and SPAM(l) behave virtually the same, and need less outer iterations than the Lanczos method to 
arrive at a given accuracy. 

Conclusions The main goal of the experiments so far was, first of all, to investigate if Full SPAM 
is competitive with the Lanczos method in terms of the number of iterations in the outer loop. 
Solving the eigenproblems for Akin the inner loop, even if this would not be done to full precision. 
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Figure 6. Lanczos versus Full SPAM and SPAM(l) for bcsstk04, bcsstkOV and bcsstklO, with approximating 
matrices from below. All three methods used the same start vector. 



makes SPAM always more expensive than Lanczos. The experiments in Figures 3, 4, 5 and 6 show 
that this is the case for different matrices A and different types of approximations Aq. Not only 
for the largest, but also for other eigenvalues of A. Secondly, we have illustrated that the use of 
one step of the Jacobi-Davidson method to approximate the eigenvalue problem for A^ in the inner 
iteration for SPAM hardly influences its behavior. We will now investigate what happens if the linear 
system in SPAM(l) is approximated with a small number of steps of the Minimal Residual method, 
MinRES [18]. 



4.5. Lanczos versus SPAM(1/): approximating the correction equation ofSPAM(I) using MinRES 

In this section we investigate the use of £ iterations of the MinRES method [18] for approximating 
the solution of the Jacobi-Davidson correction equation (33) for the eigenproblem for Ak in the 
inner iteration of SPAM(l). Each iteration of MinRES requires one matrix vector product with A^, 
which, for small values of /c, will be approximately as costly as a matrix vector product with Aq. 
The initial vector of all methods to which the comparison is applied, i.e., the eigenvector of interest 
of Aq, will still be computed in full precision. The resulting method is abbreviated by SPAM(1,^). 

Remark 4.2 

In SPAM(1,1), MinRES produces an approximation of the Jacobi-Davidson correction equation (33) 
for Akin the one-dimensional Krylov subspace with the right-hand side of that correction equation 
as start vector. Since this is the current residual from the outer iteration, the expansion is the same 
as for the Lanczos method. We will therefore not display the convergence curves of SPAM(1,1). 

In the light of the previous remark, it is reasonable to expect that SPAM(1,£) will represent 
a transition between Lanczos and SPAM(l). Thus, for reference, in the experiments to come, 
we displayed the convergence graphs of Lanczos and SPAM(l) in solid black lines without any 
additional symbols. The experiments concern four of the situations that we have already studied. 
First, we approximated the smallest eigenvalue of the reaction-diffusion problem. For the other three 
experiments we approximated the largest eigenvalue: in the second, we took the banded matrix with 
low rank approximation from below, and in the third and fourth the matrices bcsstkOV and bcsstklO 
with the respective approximations from the previous section. The results are displayed in Figure 7 
and confirm the expectations. Even for ^ = 2 and ^ = 3, SPAM(1,^) resembles SPAM(l) much more 
than it resembles Lanczos. It depends, however, very much on the actual application if the gain in 
the number of iterations is not undone by the costs of £ steps of MinRES per outer iteration with 
the matrix A^. For instance, in the banded matrix example (second picture in Figure 7), the matrix 
A itself has 322 nonzero entries and is of full rank, whereas Aq only has 33 nonzero elements 
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and is of rank 6, and especially when k is small, the action of A/^ will not be very much more 
expensive than the action of Aq. This, however, brings us to the next question that we would like 
to investigate, which is, if using in the k-th inner iteration instead of Aq all along, is going to 
make any difference, because this is what distinguishes SPAM from Jacobi-Davidson with one- step 
approximation as preconditioner. As argued in Section 2.5.4, if SPAM is going to to better, then 
probably not by very much. 




Figure 7. Lanczos and SPAM(l) compared with SPAM(1,^) for small values of i. Left: reaction-diffusion 
problem, smallest eigenvalue. Other pictures: largest eigenvalue. Middle left: banded matrix; middle right: 
bcsstkOV; right: bcsstklO. The graphs for Lanczos and SPAM(l) can also be found in Figures 4, 5, and 6. 



4.6. Comparing SPAM( with one-step preconditioned Jacobi-Davidson 

In the /c-th inner iteration of SPAM(1,^), the Jacobi-Davidson correction equation (33) for A^ is 
solved using I steps of MinRES. We will now compare this with the Jacobi-Davidson method with 
one-step approximation as preconditioner, as was described in Section 2.5.2. This means that in 
each inner iteration, the initial approximation Aq is used instead of A^. We will still apply I steps 
of MinRES to solve the corresponding correction equation and denote the resulting method by 
JD(1,£). Since one of the aims of SPAM was to save on the costs of the matrix- vector products in 
the inner iteration, we will also apply Jacobi-Davidson without preconditioning, and approximate 
the exact correction equation (29) with matrix A by performing I steps of MinRES as well, and 
denote this method by JD(^). Thus, the differences between these three methods lie in the inner 
iteration: SPAM(1,^), JD(1,^) and JD(£) all apply I steps of MinRES per inner iteration step, to a 
linear equation with matrix Ak, Aq, and A, respectively. Note that in [26], no explicit comparison 
of SPAM with Jacobi-Davidson was made, even though the methods are so closely related. 

As expected, JD(£) is the clear winner in all experiments, although the difference with JD(1,^) 
and SPAM(1,^) is not enough to automatically disqualify the latter two. Since the matrix- vector 
products in their inner iterations are, in general, considerably cheaper than in JD(^), both methods 
could be competitive. Having said this, the difference between JD(1,^) and SPAM(1,^) is quite small 
and not always in favor for SPAM(1,^), even though SPAM(1,^), much more so than JD(1,^) uses the 
information that is available from the outer iteration also in its inner iteration. As argued already in 
Section 2.5.4, this may actually be less effective than using only the best information that is available 
from the outer iteration, as the Jacobi-Davidson method does. So far, the numerical experiments are 
in favor of using Jacobi-Davidson with one-step preconditioning instead of SPAM(1,^). 
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5. CONCLUSIONS 



The experiments above illustrate mathematical aspects of SPAM as a method for approximating 
eigenvalues of a Hermitian matrix. Using approximations from below, SPAM can be seen as a 
boosted version of the Lanczos method in the sense that convergence towards the largest eigenvalues 
is stimulated. Since Lanczos itself is often used to provide a good start vector for Jacobi-Davidson, 
SPAM is a therefore a good candidate for this task, too. Since the difference between SPAM 
and Jacobi-Davidson with one step approximation is small, it may be preferred to use the latter, 
especially since the latter is even more easy to use. There does not seem to be a significant gain in 
re-using the action of A on the orthogonal complement U of the current Ritz vector u within V also 
in the inner iterations in comparison with only re-using the action of A on u, as Jacobi-Davidson 
with one step approximation does. This does not mean that the original idea of the authors [26] of 
SPAM, to save on the costs of the inner iterations of for instance Jacobi-Davidson, was incorrect. It 
may well pay off to do so, but this may be done with Jacobi-Davidson with one step approximation 
just as well. Thus, the main conclusion of this paper is that the value of SPAM probably lies in 
providing good initial approximations for the Jacobi-Davidson method. 
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